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qh. Abstract 

U ■ 

W)' A novel analytical method for f{R) modified theories without matter in Fried- 

mann-Lemaitre-Robertson- Walker spacetimes is introduced. The equation of mo- 

'sj" , tion for the scale factor in terms of cosmic time is reduced to the equation for 

J^ I the evolution of the Ricci scalar R with the Hubble parameter H. The solution 

f\j ' of equation of motion for actions of the form of power law in Ricci scalar R, is 

CN . presented with a detailed elaboration of the action quadratic in R. The reverse 

; I use of the introduced method is exemplified in finding functional forms f{R) which 

\— ^ ■ lead to specified scale factor functions. The analytical solutions are corroborated 

^vq . by numerical calculations with excellent agreement. Possible further applications 

^-H I to the phases of inflationary expansion and late-time acceleration as well as f{R) 

^ ' theories with radiation are outlined. 

c3 ! 1 Introduction 

The use of General Relativity [T] enables us to describe and explain many features 
and observational facts about the Universe we live in through enormous periods 
of time. Standard model of big-bang cosmology is a very good description of the 
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evolution of our Universe. There are of course unresolved puzzles, for instance flat- 
ness or horizon problems, which could be explained by some version of inflationary 
model OEHHE] which would provide our Universe with an early-time acceleration. 
This phase should occur before the radiation dominated epoch of cosmic evolution. 
Speaking of accelerated phases in the life of the Universe, there is also another one, 
that is a late-time acceleration following a matter dominated epoch [GlITllHllQlllOj. 
There are several observational indications of this phase of which the supernova la 
results [111 \T7[ [13] are probably the most well known example. The source of this 
late-time acceleration is frequently identified with an exotic type of matter/energy 
called dark energy and there have been many ideas and attempts in the literature 
trying to explain its origin |14[ [15] . 

In order to accommodate the fact of the late-time acceleration of the cosmic ex- 
pansion, we need to make some modifications in the standard cosmological model. 
Apart from the already mentioned possibility of the existence of dark energy, a 
prominent possibility is that the gravitational interaction is modified at (at least) 
cosmic scales. One possible way to modify General Relativity are the so called 
f{R) theories [T6l[ni[M[ia[2a[2ll[22l[2a[Ml[25l[Ml[271[28l[29l[30l[a^ where 
some function of R is used as a Lagrangian density. Starobinsky [2] used one of 
these, with a specific form given by f{R) = R + aR? to propose one of the first in- 
flationary models. f{R) theories could possibly describe the late-time acceleration 
phase, thus providing an effective model for dark energy. It is a demanding task 
to produce such a model without spoiling the successes of the standard big-bang 
cosmology while at the same time satisfying solar system and other astrophysical 
constraints. A number of such models has been proposed with various observable 
consequences and implications [321 1331 IMl 133 IMl 133 I3S] • Modified gravity theories 
have been also found useful in the resolution of the cosmological constant problem 

[32 sni mi na Ha El ss] ■ 

The principal goal of this paper is to present a novel method for analytical 
solutions of equations of motion arising in f{R) theories in FLRW cosmologies 
without matter. In theories of the f{R) type, the equation of motion is a third 
order differential equation for the scale factor a in terms of cosmic time t. The main 
virtue of the proposed method is that the said third order equation is decomposed 
into three first order ones which can be solved consecutively, i.e. one by one and 
not simultaneously as a system of first order equations. 

The organization of the paper is the following. The first section is the introduc- 
tion. In the second section the formalism of f{R) theories is briefiy summarized 
and analytical approaches and solutions of f{R) theories in the existing literature 
are discussed. The third section contains the description of the novel method. In 
the fourth section the introduced method is applied to theories f{R) ~ i?" with a 
special emphasis on actions quadratic in R. The fifth section describes the reverse 
application of the method: finding the functional form of f{R) which results in a 
specified scale factor function, a{t). In the sixth section we present a numerical al- 
gorithm used for the verification of analytical solutions and in the seventh section 
we present the comparison of analytical and numerical solutions. In the eighth 
section we describe how to extend our method to f{R) theories with radiation. 
The paper closes with the conclusions section. 



2 The f{R) modified gravity theories in FLRW 
cosmologies 

In f{R) theories the gravitational part of the action is 

S^r.. = -^jd^x./^gf{R). (1) 

The variation of the entire action S = /Sgrav + 5'matter , where 5'matter represents the 
matter part of the action, with respect to the metric g^^ yields equations of motion 

f{R)R^u - \f{R)g^u + (V^V, - gf,,n)f'{R) = -SttGT^, . (2) 

In the FLRW metric 

ds^ = dt" - a{tf (j^^ + r^dnA , (3) 

the 00 component of ([2]), for the spatially flat 3D space, then reads 

3f'iR)H^ - ^{Rf{R) - f{R)) + SHRf'iR) = SnGp . (4) 

Despite widespread numerical methods and computational power nowadays 
available for finding numerical solutions of the equation of motion Q, the most 
complete insight into the functioning of any theory is achieved when analyti- 
cal solutions can be obtained, f{R) theories not being exception to this rule. 
There have been many attempts of finding analytical solutions for f{R) theories 
[l6l[I71[ig[ia[20l[2ll[22l[2a[2a[23[2S[271[Ml[l9l[3a particular, the cos- 

mological power law solutions for the scale factor dynamics in universes with and 
without matter were found in [IBIHT] (see also [l8j) for f{R) ~ -R". The existence 
of Godel, Einstein and de Sitter universes in modified gravity theories was studied 
in [49j whereas analytical solutions for spherically symmetric systems were found 
in |50t l51j. Methods of finding analytical solutions of f{R) theories using Noether 
symmetries were presented in |52| . 

A systematic approach to solving f{R) equations of motion was adopted in [53] • 
Using a specific scalar field reformulation of f{R) theories, Clifton succeeded in 
finding analytical solutions for R^'^ theories in vacuum spatially flat FLRW spaces. 
Furthermore, analytical solutions for spatially curved vacuum FLRW spaces and 
some specific solutions with perfect fluids were also found in j53j . The method 
used in [53| requires a number of quite ingenious, but intricate substitutions and 
changes of variable. 

In this paper, we present a general approach to finding analytical solutions 
in f{R) theories. The proposed method reduces the third order equation into 
three first order ones which can be solved one after the another, i.e. do not have 
to be solved simultaneously as a system of first order equations. We also present 
solutions of i?" theories using our approach as a demonstration of the method. The 
selected models of the type ~ R" were primarily chosen for illustrative purposes, 
especially because they have so far attracted a lot of attention of researchers in the 



field of modified gravity. Tliese models have been so far mostly ruled out by the 
observations. In the light of results of [53], our results for i?° theories cannot be 
considered novel, but their alternative derivation in this paper, simpler and more 
straightforward in our opinion that that in |53j . is new. The greatest virtue of the 
introduced method is its simplicity and generality. 

In the remainder of the paper we restrict our considerations to spatially flat 
spaces {k = 0) without matter, i.e. to vacuum solutions in flat 3D spaces. Only in 
section 8 we discuss the addition of radiation in f{R) theories. 

3 The method 

The dynamical equation describing the expansion of the spatially flat universe 
without matter (or when the contribution of matter is negligible) is 

3FiR)H^ = ^{F{R)R - f{R)) - 3HF{R) , (5) 

where F{R) = f'{R), prime denotes the differentiation with respect to R and dot 
represents the differentiation with respect to cosmic time. The crucial step of the 
method presented in this paper stems from the structure of the Ricci scalar in the 
RW metric: 

R = 6H + 12H^ . (6) 

Namely, in intervals of cosmic time in which H{t) is a monotonic function of t, it 
is possible to make the change of variables from t to H. In particular. 
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(7) 

If H{t) is not a monotonic function of cosmic time t, the time domain should be 
divided into intervals where H{t) is a monotonic function and the transformation 
([7]) should be performed in each interval separately. Inserting d?]) into ^ yields 

3H (I - 2H'^ f"(R)^ = lif'iR)R - fiR)) - ^H'f'iR) . (8) 

The representation of the starting third-order differential equation ([5]) in the form 
given above opens a way for a program of obtaining analytical solutions by solving 
three first-order differential equations one after the other. In the first step, the 
solution of ([5]) results in the function R{H). Once the Ricci scalar R is available as 
the function of the Hubble expansion rate H, solving ([6|) results in H{t). Finally, 
knowing H{t), by the definition of the Hubble expansion rate, leads immediately 
to a{t). 

Which of the three steps described above can be carried out analytically de- 
pends on the very form of the function f{R). It should be stressed that the crucial 
step in ([7|) is not possible for geometries with k ^ 0. The remainder of this pa- 
per is devoted to the study of cases where at least one of the steps can be done 
analytically. 



4 The applications 

4.1 Action quadratic in R 

The action quadratic in Ricci scalar is of physical interest for several reasons. 
Prom a purely formal point of view, it is the power-law term of the lowest positive 
integer power which yields a nontrivial f{R) modified gravity theory. The terms 
quadratic in R appear as a part of vacuum action in quantum field theory in curved 
space-time, see e.g. [5l]. Finally, R^ terms play a crucial role in one of the most 
prominent models of inflation, the Starobinsky model [2J (i?^ term is actually a 
truncation of the Starobinsky model f{R) = R + AR^ and it becomes important 
in the early universe when the curvature is large) and in some recent analyses of 
the radiation epoch of the cosmic evolution [55]. 

Inserting the function f{R) = AR^, where j4 is a real constant, into ([8]) gives 
the equation 
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H{R-12H')— = -R{R-12H'), (9) 



or more compactly 
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{R-UH'){H—--R]=0. (10) 



As long as H is evolving in time (and R — 12H ^ 0) the equation to be solved is 
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The solution for R{H) is then 
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R{H) = ^0 ( ;^ ) , (12) 



where R{Hq) = Rq is the initial condition. 
Eq. ([6]) now acquires the form 

dH Ro f H 



1/2 
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with the following closed-form solution: 



6^2/3 ^2/3 + ^i/3{H/Ho)y^ + {H/Ho) 
+ ^arctan^^^^^^^l^^^±^ 

1 ^ ix'/' - If 1 , 2 + x'/' ,,,, 

+ ^^ ^^ x^'^ + x'/' + 1 " ^^Wl V^^ ' ^ ^ 

where x = Ro/{^'^Hq). This solution implicitly defines H in terms of t. From the 
definition of the Hubble parameter, H = a/a, one readily obtains the equation for 
the evolution of a with H. Namely, from (|13p it follows 

da HdH H dH 
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A straightforward integration of this equation yields 



a = ao 
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whereas an inversion of this expression yields 



H = Hn 



X - (x - 1) 
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(16) 
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(17) 



The expression (J17p was already obtained in [56j using a different method. This 
solution can be inserted into ()14p to obtain the implicit relation connecting t and 
o. Finally, from the solutions ([H]) and (fT7|) one can easily read out the asymptotic 
behavior: for t — t- oo, the scale factor is unbounded, a — )■ oo, whereas the Hubble 
parameter saturates, H — )■ x^ Hq. 



4.2 Separability oi R - 12H^ term 

The effectiveness of the method for the action quadratic in R is closely connected 
to a very simple form of equation for i? as a function of H, Eq. (|11|) . This simple 
form is a consequence of the factorization of the R — 12H^ term or the right-hand 
side of ([8]), as explicitly demonstrated in pO|) . An interesting question is which 
other, if any, functional forms f{R) also allow the factorization of the R — 12H^ 
term on the right-hand side of dH). Indeed, the factorization 
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U'{R)R - f{R)) - ^H^f'iR) = t{H)X{R){R - 12H^ 



(18) 



where t{H) and X{H) are arbitrary functions of H and R, respectively, results 
in a separable equation once the R — 12H^ term is factored out. Term by term 
comparison of right and left-hand side in (|18|) gives t{H) = const which can be 
absorbed into X{R) so that we can take t{H) = 1. Further we obtain 



and 



X{R) = -^f'{R) 



-(f'iR)R-fiR)) = XiR)R. 



(19) 



(20) 



Combination of (J19p and (j20p gives f{R) = AR^. Therefore, the action quadratic 
in R is the only choice for f{R) which allows factorization as described in ()18p . 



4.3 f{R) = AR'' 

The general power law form f{R) = AR°' is a natural extension of the quadratic 
action. This choice of action has recently raised a lot of attention in the literature 
as adding such a power law term to the R term in the modified gravity action 
can lead to accelerated late-time expansion of the universe such as for negative 
values of a. The action f{R) = AR" is chosen primarily for illustration purposes. 



It should be stressed that it is strongly constrained [IHl [50| and disfavored as a 
model of late-time universe by light-deflection and planetary period observations 
|57| and the requirement of avoiding the Dolgov-Kawasaki instability [58]. 
For this choice of f{R) Eq. ^ acquires the form 

(77? 

a(a - 1)(R - 12H'^)H = (a - l)R^ - GaH'^R . (21) 

dH 

Using the substitution R = ^H^, the equation above can be transformed into a 

separable form 

dH _ a e-12 , 

_ 6a(4a - 5) 
^ -(«_!)(!_ 2a)- ^^"^^ 
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Writing 



12 _Ci ^ C2 _ (Ci + C2)g + gi/3 ^24) 



?(^ + /3) e e + /3 e(e + /3) 

we can see that 

Ci + C2 = 1 , 

^'""7' (25) 

Equation ()22p can be readily integrated to obtain the solution in closed form: 

In general, Eq. ()26p can be considered as a parametric solution of the R{H) 
relation. Namely, (j26p yields H{(^), whereas R{£,) is obtained directly from its 
definition as R{0 = ^i?(0^- 

Let us further consider some values of a for which (|26p can be inverted to obtain 
the explicit relation R = R{H). Defining 

,..«£%±^, (27) 

([26]) can be written as 

C^' (^ + /5)^' = KH^ . (28) 

Let us examine the following combinations of Ci and C2 (and therefore values of 
a): 



1. Ci = — 2C2 will give C2 = — 1 and Ci = 2 and the equation ()28p becomes 

^'(? + /3)-' = i^^'^, (29) 

with solutions 



(1^2 = -KH^ {l±\/l-24K-^H-^ir' ] . (30) 



This leads to solutions for R 



Ri,2 = ^H^{l±\/l-24K-^H-^] , (31) 



where 



"1,2 = — ^ — • (32) 

2. C2 = — 2(7i will give Ci = —1 and C2 = 2 so this time ()28p becomes 

^-i(e + /3)2 = i^F^, (33) 

leading to solutions 



2 

so that 



^12 = - I -24 + KH^ ± J ( KH"^ - 48 ) KH"^ ) , (34) 



i?i_2 = -i?^ ( -24 + KH^ ±\l( KH^^ - 48 ) KH^ 



\ 1-2q \ 

j K//^^ , (35) 



where 



ai,2 = Yq • ^^^' 



3. Another case which would lead to quadratic equation in ^ would be Ci = 
C2 = ^ = — ^ where we would have /3 = — 24 which would lead to a complex 

values of a: 

7±^/^T5 ,„_, 

"1,2 = ^ • (37) 

Finally, the inspection of ()26p and the definitions of Ci, C2 and /?, reveal 
the following characteristic intervals and point values of a: (00, 0), [0,1/2), 1/2, 
(1/2,1), [1,5/4), 5/4, (5/4,2) and [2, 00). For points of q = 1/2,5/4 we solve ([21]) 
directly by inserting concrete values of a. This results in 
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and 

H = Hoe''^^ ^-j , (39) 

for a = 1/2, 5/4 respectively, with R = ^H^. 

5 Functional form of f(R) from the known so- 
lutions for the scale factor 

The preceding section was dedicated to finding analytical solutions for the known 
functional form f{R)- The ordering of this approach can also be reversed. Namely, 
for a known scale factor dependence on time, a{t), one may ask which functional 
forms f{R) yield such a{t) functions as solutions. We consider the reconstruction 
in the mathematical sense, i.e. obtaining the function f{R) which without matter 
leads to a known function a{t). The procedure is the following: from a known 
a{t) calculate a{t), then H(t) and H{t). Using ([G]), the expression for R{t) follows 
immediately. If from expressions for H{t) and R{t) it is possible to eliminate t 
and obtain relation H^ = k{R), the relation ([8]) becomes a first order differential 
equation for the function f{R). For related work on the reconstruction methods 
in f{R) gravity see [591 [601 ED [62l [Ml El] • 

As an illustration of the procedure explained above, we determine functional 
forms f{R) which lead to power law expansion, a{t) ~ t^ and singular future 
behavior a{t) ~ 1/{T - t)"". 

5.1 Power law expansion 

For the power law expansion 

a{t) = Dt^ , (40) 

where D and /3 are constants, we have H{t) = j3/t and H{t) = —fi/t^. The 
expression for the Ricci scalar is R{t) = 6/3(2/3 — l)/t^ which is readily displayed 
as i? = ^ g ' H^ = -H^. Eq. ([8|) now acquires the form 

(1 - l2^)R^f"{R) - i(l - 67)ii/'(i?) + \f{R) = . (41) 

This equation has solutions of the form f{R) ~ R^ where A is the solution of the 
equation 

2A2 + (/3-3)A + 1-2/3 = 0, (42) 

in accordance with the results of |i46f i47j. Direct inspection of this equation shows 
that A = 2 is not its solution for any value of /3. Solving ([^2]) follows 

1 



Ai,2 = |(3-/3±v'/5' + 10/3 + 1). (43) 

In an expanding universe /3 > and Ai 7^ A2. This result shows that the functional 
form f{R) that leads to the power law expansion (|40p is a linear combination of 
R^^ and R^^ . See also [651 for a similar result with matter. 



5.2 Future singularities 

As some mechanisms for the explanation of the late-time accelerated expansion of 
the universe lead to singularities in finite time in the future |29] (phantom energy 
with the constant equation of state parameter being one such example), study of 
such expansion scenarios is of interest in modified gravity theories as well. For a 
scale factor evolution with a future singularity at time T > Iq 

^ (44) 



(T - t)"^ ' 

where A and m are positive constants, the Hubble parameter and its derivative 
have the form 

T-t {T-ty ^ ' 

The Ricci scalar is again quadratically dependent on H: 

R^'J^ail^HK (46) 

m 

Using this result ([S]) becomes 

2R^f"{R) - (m + l)Rf{R) + (2m + l)/(-R) = . (47) 

The solutions of this equation are in the R form where the exponents A are 

TTi + 3 lb y/m? — lOm + 1 
Ai,2 = 1 . (48) 

For vfi < mi = 5 — 2\/6 and m > m2 = 5 + 2\/6, both solutions for A are real and 
the solution for f{R) is 

f{R) = KiR^' + K2R^^ , (49) 

where K12 are real constants. For m = mi^2 the solution for f{R) is a linear 
combination of i?(™i>2+3)/4 ^^^ i?(™'i.2+3)/4 j]^^_ Finally, for mi < m < m2 the 
solutions for A are complex conjugate, Ai^2 = Ar it A/, with Xr = ^^^ and Xr = 
-^^J\vn? — 10m + 1|. The solution for fiR) leading to ()44p is 

/(i?) = K^R^R cos(A/ In R) + K^R^^ sin(A7 In R) , (50) 

where ^^3,4 are real constants. The obtained results are consistent with |29j. 

The examples discussed in this section illustrate one more additional advan- 
tage of the analytical method introduced in this paper. It allows us to find some 
/(i?) actions leading to predefined function a(t). Thus, it could potentially add 
to the existing literature on the reconstruction of modified gravity models which 
mimic known expansion epochs such as inflation, matter dominated or radiation 
dominated epoch of the expansion of the universe. 
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6 Numerical solutions 

When we have knowledge of the f{R) function, we want to find the evolution of 
the scale factor and other quantities (Hubble parameter, curvature...). As we 
have seen in the above sections, unlike standard Friedmann equations, analogous 
equations for arbitrary f{R) are very complicated and hard to solve. When the 
matter is absent the equation requiring solution is 



9 / • • • 

a \ 1 ,^,„^„ „,„XN a dF a da 




mR) [-)=-, (F(R)R - m)) -is--^---,-2{-:] I . (51) 

It can be rewritten in a more suitable form for numerical integration as a system 
of three differential equations of the first order 

d = b, (52) 

b = a = c, (53) 



For a general class of /(i?) functions 

f{R)=PR + a{R-Ro)'', (55) 

the parameters a, f3, Rq and n can be chosen in such a way to recover some theories 
where the solutions are known. The only restriction is not to go to standard GR 
plus cosmological constant limit, because then we would have n = or n = 1 or 
a = because this would lead to divergences in (j54p . Parameter Rq serves to avoid 
possible divergences in the n < case. 

For f{R) chosen as in ([55]) equation (j5l|) becomes 



\n-2 



18ban{n-l) {R - Rq 
- 3{/3 + aniR-Ro)''-^)^)--(c-—) (56) 

(Jj / (Ji \ (Jj / 

The output is given in terms of a{t), b{t), c{t) functions so that any other 
quantity we need can be calculated, such as R{t), H{t), or q{t). This allows 
analysis of functional dependencies, such as R{H). 

We need to make our variables dimensionless. To accomplish that take r = 
HQ{t — to) leading to dr = H^dt. Also let a = a^x where ag is the scale parameter 
today. Then 

d = aoHl^. (57) 
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Any other quantity is scaled accordingly, so for example 



TT TJ 1 dx 



R = QHl 



X dr 
ldx\'^ I ld^2 



_(^f) +^iij (58) 

and dimensionless quantities can be constructed 

H R , . 

h = -Er^ r = T72- (59) 

-no ^0 

The graphs representing numerical solutions generated by this program along 
with analytical solutions for several different choices of f[R) = -R" are discussed 
in the following section. 

7 Results and discussion 

In this section we verify our analytical solutions for i?° theories obtained in section 
4 by comparing them with numerical solutions obtained using algorithm laid out 
in section 6. The results of the comparison for all characteristic point values and 
intervals for a are presented in Figures [TMl The line in the figures represents a 
parametric plot of R{^)-H{S,) which is obtained using (|26|) and R = S,H^ while 
the dots represent the output of the numerical procedure. We reiterate that it is 
important to apply our analytical method only in intervals of cosmic time in which 
H{t) is a monotonic function. 

A distinctive feature of some of the plots presented in Figures [T][6] is the existence 
of extrema of functions -ff (^) and/or -R(^). In order to analise the extrema of R{^)- 
H{S,) functions we write 



H = AC (e + /3)^' (60) 



and 



R = H^^ = Bf^^+\^ + l3f''\ (61) 

where 71 = — 2 ^^~ ^ and 72 = (4a-t)(i-2a) • -^^o™ these two equations we have 

^ - f - :^f7.+i^. + ,^. j27i + l)(^ + /?) + 2^72 . . 

dH-f-A^ ^^ + ^^ 7i(e + /3) + e72 • ^'^ 

Looking at this expression we can see that R has an extremum for ^ = if 
71 + 1 > (equivalently a £ (—00, 5/4) |J (3/2, 00) ) and the same goes for H 
if 71 + 1 < (for a G (5/4,3/2)). Also at ^ = -/3 for 72 > (equivalently for 
a € (—00, 1/2) IJ(5/4, 2)) R has an extremum while H has an extremum for 72 < 
(for a G (1/2,5/4) U (2, 00)). Also the analysis of possible extrema coming from 
the numerator and denominator of the fraction in Eq. (|62|) one can see that the 
condition for R to have an extremum is f = -^ and for the extremum in H 
condition is .^ = 12. 
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Figure 1; Functional dependence of dimensionless quantities r(h) for f(R) = B? 





Figure 2; r{h) dependence for f{R) = R ^ and f{R) = R 
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Figure 3: r{h) dependence for f(R) = R°'^^ and f(R) = R° 
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8 f(R) theories with radiation 

The method of solving f{R) theories can be extended beyond vacuum spaces. To 
do that, it is necessary to go beyond the 00 component equation of motion. The 
trace of the equation of motion ([2]) yields a relation 



Rf\R) - 2f{R) - 3nf{R) = -SttGT , 



(63) 



where T = Tf^ is the trace of the energy-momentum tensor of matter. In the FLRW 
cosmologies the d'Alambertian operator can be represented as 



--w-< 



(64) 
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Figure 4: r{h) dependence for f{R) = B?-'^^ and f{R) = R^-^ 
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Figure 5: r{h) dependence for f{R) = R^/'^ and f(R) = R^-^ 
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Figure 6: r{h) dependence for f{R) = R'^ and f{R) = R* 

In combination with the expression ^ = ("f — ^-^^)d7r' ^^ ^^ evident that in FLRW 
the d'Alambertian operator can be expressed in terms of derivatives with respect 
to H. On the other hand, if the matter is dominated by radiation, the trace of 
energy-momentum tensor T = pr — 3pr vanishes. The equation of motion for f{R) 
theory with radiation can therefore be represented as 



Rf'{R)-2f{R)-3(^^-2H''^^ 
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df'jR) 
dH 



-9H{-- 2H' 



df'jR) 
dH 
(65) 

As in the vacuum case, the problem is reduced to three independent differential 
equations which can be solved consecutively. The case with radiation, however, 
differs from the vacuum case in the fact that the differential equation for R{H) is 
of second order. 



0. 



14 



Apart from a direct application on f{R) theories with radiation, Eq. ()65p 
is also applicable in the case of no matter, since T = in that case too. The 
method represented by Eq. ([SS]) may, therefore, also be applied in the study of the 
transition from the vacuum case to the radiation case or vice versa, should this 
case appear in some f{R) model. The transition from the epoch of inflation to the 
epoch of radiation might be the first place where the potential of ()65p might be 
tested. 

9 Conclusions 

We have presented a novel analytical approach to solving vacuum equations of 
motion for f{R) theories in FLRW spaces. A key advantage of the method is 
the decomposition of the third-order differential equation into three first-order 
ones which can be solved consecutively. Compared to methods used to obtain 
some exact solutions in the literature [S3], the approach introduced in this paper 
is simple and much more universal. The introduced method also provides an 
approach to finding functional forms f{R) which produce a given dependence of 
the scale factor on cosmic time. This fact paves the way for the use of the method in 
building f{R) models which closely mimic some important phases in the evolution 
of the universe and which could hopefully be useful in tackling problems of dark 
matter and dark energy. We also outline the application of the method to f{R) 
theories with radiation. As vacuum solutions of f{R) theories in FLRW spaces are 
important in the study of early (inflation) and late-time universe (once the matter 
density is sufficiently diluted by the expansion), we hope that the universality and 
simplicity of the proposed method will make it useful in future applications of f{R) 
theories in cosmology. 
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